Magnetism and pairing of two-dimensional trapped fermions 
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The emergence of local phases in a trapped two-component Fermi gas in an optical lattice is 
studied using quantum Monte Carlo simulations. We treat temperatures that are comparable or 
lower than those presently achievable in experiments and large enough systems that both magnetic 
and paired phases can be detected by inspection of the behavior of suitable short-range correlations. 
We use the latter to suggest the interaction strength and temperature range at which experimental 
observation of incipient magnetism and d-wave pairing are more likely and evaluate the relation 
between entropy and temperature in two-dimensional confined fermionic systems. 

PACS numbers: 



Strong electron correlations in solids are believed to 
be at the origin of a remarkable host of phenomena that 
range from anomalous insulating states and magnetism 
to high temperature superconductivity [1]. The Hub- 
bard Hamiltonian [2] , a model which provides a simplified 
picture of the band structure and electron interactions, 
is thought to contain the necessary ingredients to de- 
scribe such spectacular diversity. In the strongly interact- 
ing limit the half-filled model describes a Mott insulator 
and, on a two-dimensional (2D) square lattice, quantum 
Monte Carlo simulations [3, 4] have convincingly estab- 
lished the existence of antiferromagnetism at T = 0. This 
conjunction of insulating behavior and long-range anti- 
ferromagnetic order accurately describes the low temper- 
ature physics of the undoped parent compounds of the 
cuprate superconductors. 

Despite intensive analytical and computational work, 
what happens as one dopes the antiferromagnet has re- 
mained controversial and many important questions re- 
main unanswered [5-7] . A promising new route to study 
the model has recently emerged in the form of experi- 
ments with fermionic gases in optical lattices. The ap- 
peal of these experiments lies in the high degree of con- 
trol over the different parameters of the system and the 
ensuing possibility of a close realization of a "pure" Hub- 
bard Hamiltonian [8] , free of the complexity characteriz- 
ing solid state systems. Ground breaking achievements 
include loading an ideal quantum degenerate Fermi gas 
in a three-dimensional (3D) lattice [9] and the realization 
of the Mott metal-insulator transition in 3D [10, 11]. 

Optical lattice experiments require the presence of a 
confining potential, usually harmonic, and are accord- 
ingly modeled by the Hamiltonian 

^ = E (4'^- + H-C.) + + Un^tn,^) , (1) 

with Vi = Vrf. The curvature V causes the density to 
vary across the lattice, resulting in a situation at odds 
with the homogeneous system that one would like to emu- 



late. This is a potential problem because numerical stud- 
ies on related models [12] have indicated the existence 
of many competing phases with vastly different physical 
properties. Understanding the extent to which the deli- 
cate balance between these phases is affected by the pres- 
ence of the trap is, therefore, of importance for assessing 
whether fermionic gases can be taken as accurate simu- 
lators of homogeneous models. Thanks to improvements 
in quantum Monte Carlo (QMC) codes [4, 13-16], this 
and other issues can now be quantitatively investigated 
at temperatures that are comparable or below those of 
the latest optical lattice experiments [10, 11]. 

Here we report results from state of the art determi- 
nant QMC (DQMC) simulations [17] of the trapped 2D 
Fermi-Hubbard Hamiltonian. We found clear signatures 
of magnetic and pairing correlations at temperatures of 
the order of the magnetic scale J = At'^/U, a perhaps 
surprising observation considering that the correspond- 
ing temperature range in strongly correlated solid state 
systems is well above Tc- We note, however, that the 
properties that we compute are local in character, and 
that recent advances in optical lattice experiments [18- 
20] have made possible the imaging of individual sites and 
short range correlations around them. Hence, our results 
suggest that the purity of optical lattices and the novel 
probes used in these experiments could allow the obser- 
vation of local pairing signatures at higher temperatures 
than possible in a condensed matter environment. 

We analyze the properties of the trapped system in 
terms of the spatial dependence of the local density 
n{i) = (uj), the density fluctuations 6{i) = (nf) — (n^)^, 
and the spin correlations Cxy{i) = 4(5'^'^^^ y)^i'^' ^^^^^^^S 
is analogously discussed by defining a local d-wave pair 
creation operator A^^ 

f (2) 

A"*^ = -(A^^^ + At . .-A^^.-A^ o 
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FIG. 1: (a) Local density n{i), (b) density fluctuations S{i), 
(c) local staggered magnetization M{i), and (d) d-wave pair- 
ing are shown for U = 6t, T = 0.31t, /io = 3.0t, and 
V = 0.04t. A Mott insulating domain is emerging in the den- 
sity profile of panel (a), in the form of a half-filled ring 6-10 
lattice spacings from the trap center. The density fluctua- 
tions are minimized in this region. The staggered magneti- 
zation and the d-wave pairing, however, show a pronounced 
maximum in the Mott domain. 




FIG. 2: Temperature dependence of the (a) density, (b) 
density fluctuations, (c) nearest-neighbor spin correlations 
Cio(i), and (d) next-nearest-neighbor d-wave pairing -Pii(i) 
for U = 6t. Lines are results within the local density ap- 
proximation. The insets of (a) and (b) show the temperature 
dependence of the entropy and the double occupancy normal- 
ized to the number of particles, respectively, while the insets 
of (c) and (d) are the averages over the lattice of the local 
staggered magnetization and d-wave pairing. 



and a corresponding correlation function Pxy{i) ~ 
(A^^(^ )■ To isolate the effect of the pairing ver- 
tex [21], we focus on the connected correlation function 
in which suitable products of single particle propagators 
are subtracted from the expression in Eq. (2). 

In Fig. 1, we examine the spatial dependence of dif- 
ferent properties at [/ = 6t, F = 0.04 < and T = 0.31 1 
for a system of 560 fcrmions. In this parameter regime, 
the average sign in DQMC is 0.3 and decreases exponen- 
tially as T is lowered. T = 0.31 1 is therefore very close 
to the bound that the sign problem [22, 23] imposes on 
the lowest achievable temperature. Despite this restric- 
tion, the emergence of a Mott insulator is clearly visible 
and depicted in Fig. 1(a), which shows a density plateau 
region of commensurate filling, n{i) ~ 1, and in Fig. 1(b) 
through the formation of a pronounced minimum in the 
density fluctuations. This domain is also distinguished 
by enhanced antiferromagnetism [Fig. 1(c)] and d-wavc 
pairing [Fig. 1(d)]. These last two properties are repre- 
sented by the local order parameters 

M{^) = J2i-ir^^^^^+^C,A^), X(*) = E ^-.^(*)' 

x,y x,y 

which are expected to diverge if long-range order develops 
around site i. 

Given that the temperature T = 0.31 1 is of the order 
of the antiferromagnetic exchange J, the sharp signal in 
M{i) is clearly due to the formation (and partial order- 
ing) of local moments in the Mott domain. The appear- 
ance of the peak in x(i) in the same region is however sur- 
prising - one would expect the peak to occur away from 



the insulating phase. This shows that M(i) or %(«) can be 
accurate indicators of the appearance of local phases only 
at temperatures sufficiently low that the order parameter 
is dominated by long-range contributions. As Fig. 1(d) 
demonstrates, this is not the case in our simulations nor 
in current experiments. However we shall argue that the 
temperature and interaction dependence of spin and pair- 
ing correlations, when examined at distances which are 
sensitive to the appropriate energy scales, can provide 
compelling evidence of incipient order [24] . 

Using this argument, we proceed to determine the tem- 
perature and entropy scale that experimentalists need to 
achieve in trapped 2D lattices in order to observe the on- 
set of antiferromagnetic and pairing correlations. Note 
that a related analysis has been recently carried out for 
the entropy scale of the half-filled homogeneous model 
[25]. 

As the temperature is lowered from T = 1.14 1 to 
T = 0.31 1 {U = 6<), the density distribution [Fig. 2(a)] 
is slightly compressed, the entropy per particle S [shown 
in Fig. 2(a) inset] decreases monotonically, and there is 
a marked increase in the total double occupancy in the 
system D [see inset in Fig. 2(b)], which is generated by 
the increase of the double occupancy at the trap cen- 
ter. We determine the entropy, S = /3{{E) — F), by 
first computing the Helmholtz free energy using constant- 
temperature coupling-constant integration over an iden- 
tical set of traps containing approximately 560 fermions 
and increasing U [24]. The Mott insulating region n{i) ~ 
1 is characterized by a deepening of the minimum in the 
density fluctuations S{i) [Fig. 2(b)] and by a three-fold 
increase in the nearest-neighbor (n.n.) spin-spin correla- 
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FIG. 3: Same quantities of Fig. 2 but as a function of 
interaction strength U/t at constant T = 0.50t. The in- 
set of panel (a) shows the entropy for temperatures T/t — 
0.67, 0.57, 0.50, 0.44, 0.40, 0.36, 0.33, and 0.31 (top to bottom). 



tion Cio(«) [Fig. 2(c)]. As the latter may be the easiest 
way to observe the onset of antiferromagnetic correlations 
in experiments, Fig. 2(c) makes clear that, although a 
weak signal may be visible at the highest T, a clear max- 
imum in the Mott domain will only be apparent when 
T <t. 

When analyzing pairing correlations it is important to 
realize that Poo(*) a-nd Pio{i) contain large contributions 
from Cio(i) and, consequently, carry little information on 
the actual pairing amplitude in the vicinity of a given site. 
This is analogous to the more familiar statement that lo- 
cal moment formation cannot be taken as indication of 
short range magnetic order. Pii{i) [Fig. 2(d)] is there- 
fore the shortest-distance correlation that can be used to 
gain insights on the development of off-diagonal order; it 
is characterized by two maxima around the Mott region, 
sharply peaked at densities {n = 0.8 and n = 1.2) that lie 
in the center of the superconducting dome of the cupratc 
phase diagrams [26] . The analogous correlation for pairs 
in an extended s-wave state (see [24]) is also significantly 
enhanced in a broad ring around n = 0.6. Although the 
local character of Pn and the high temperature do not al- 
low a characterization of the dominant low-temperature 
pairing symmetry, our results suggest that d-wave and 
extended s-wave pairing may be detected as leading pair- 
ing symmetries in spatially distinct regions, in agreement 
with existing ground state calculations on related homo- 
geneous models [27, 28]. 

Another important piece of information that will help 
experimentalists observe antiferromagnetism and pairing 
is the optimal value of the ratio U /t at which those corre- 
lations are expected to be maximal. In order to address 
this issue, we first investigate the interaction dependence 
of the same physical properties represented in Fig. 2 but 
at constant T = 0.5 1 (Fig. 3), and then show that our 
conclusions are unaltered if it is S to be held constant. 



At intermediate coupling {U = 4t), the smallness of 
the gap, the curvature of the trapping potential, and 
thermal fluctuations make it impossible for a flattening in 
the density profile to appear. Note, however, that signals 
of incipient insulating behavior are clear in the structure 
of density fluctuations [Fig. 3(b)]. At large interaction 
strength {U = 10t,12t), the existence of a Mott insu- 
lating domain can be directly inferred by inspection of 
the density profile: the large central plateau comprising 
about 300 sites is indicative of a phase with vanishing 
compressibility. Nearest-neighbor spin correlations ac- 
quire their largest value at U ~ 8t. While the double 
occupancy is monotonically suppressed as U increases, 
the staggered magnetization shows a similar behavior to 
Cio(j), reaching a maximum when U is of the order of 
the bandwidth. It is therefore clear that the optimal in- 
teraction strength must be near the bandwidth and that 
future experiments should focus their attention on this 
regime. Similarly Pii{i) exhibits a sharp maximum at 
Ti{i) = 0.8 for [/ ~ 8< but contrary to what is observed 
for Cio(j), there is no weakening when going from U ~ 8t 
to U ~ 12t and no dependence of the "optimal" density 
for pairing on U . The latter is reminiscent of the prop- 
erties of cuprate superconductors where optimal doping 
is constant across the entire family of these materials. 

One might wonder what would happen in the exper- 
imentally relevant case of a (quasi-)adiabatic evolution, 
i.e., a system at constant entropy, as a function oi U. In 
the inset of Fig. 3(a), we show the evolution of the en- 
tropy per particle as a function of the interaction strength 
for various fixed temperatures. As U is increased at con- 
stant S, we find that the system cools down and goes, 
for instance, from T = OM at = 2t to T = 0.3U at 
U = lot, while S is held fixed at 0.5. Since this cooling is 
mild beyond U = 8t, isothermal and adiabatic evolution 
in this parameter regime are essentially equivalent. For 
U < 8i, an adiabatic increase of U is accompanied by pro- 
gressively lower temperatures and this, in turn, implies 
that U ~ 8t remains the optimal interaction strength to 
observe strong-correlation phenomena. 

We finally touch on the extent to which properties in 
a trap are a faithful representation of those in a homoge- 
neous system. To this aim we use DQMC simulations on 
homogeneous 8x8 clusters to estimate local properties 
as a function of chemical potential and plot the corre- 
sponding local density approximation (LDA) results in 
Figs. 2 and 3 as continuous lines. The agreement with 
the ab-initio results is excellent with the sole exception 
of the half-filled region at T = 0.31 i [Fig. 2(c)]. This 
discrepancy is a characteristic failure of the LDA when 
applied to the interface between coexisting phases [29]. 

A more dramatic failure of the LDA is expected 
when examining longer range correlations. Since in the 
fermionic Hubbard model the development of a certain 
type of order corresponds with a rather narrow density 
interval, different phases in a trap appear with a narrow- 
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FIG. 4: Spin correlations Ca;y(i) [for fixed site i and variable 
separation {x, y)] along the path illustrated in inset (a) around 
the trapped center. Inset (b) is an enlarged view of the longer 
correlations. Cxy{i) exhibits the alternating sign characteris- 
tic of antiferromagnetism. The correlations decrease rapidly, 
but nevertheless extend out to several lattice spacings. The 
spin correlations are also shown in the LDA approximation 
for d = 1 and d = 2. 



netic and pairing correlations. Tuning experiments to be 
in the regime where the onsite interactions are of the or- 
der the bandwidth {U ~ St) provides the sharpest signal 
of the many-body effects. Our computation of the en- 
tropy S indicates that adiabatic cooling occurs in the 2D 
Hubbard Hamiltonian with a position dependent (confin- 
ing) potential as the interaction strength U is increased 
via, e.g., a Feshbach resonance. This allows our con- 
clusions concerning the observability of spin and pairing 
order to be relevant to experiments at constant entropy. 
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OCI-0904597 and through TeraGrid resources provided 
by NIGS under grant number TG-DMR100007. M.R. 
and R.T.S. thank the Aspen Genter for Physics for hos- 
pitality. G.N.V. acknowledges NSF dissertation enhance- 
ment grant No. 0803230. 



ring-like topology. If one disregards interfacial effects, 
these regions can be thought as one-dimensional homo- 
geneous strips where the long range character of corre- 
lations is determined by the strip width. An example 
of such a situation is given by the development of long- 
range magnetic order in the half-filled annulus found in 
our simulations a,t U ~ Qt. In particular, Fig. 4 shows 
that 2D correlations, computed on a half-filled homoge- 
neous cluster at the same temperature, represent a good 
approximation to the correlations in a trap only at short 
distances, overestimating the correct long-range value. 
Analogous results for a ID chain fail at all distances and, 
in particular, decay too quickly at large separation. This 
behavior is generic and expected for the supcrfluid re- 
gions as well. When the system size increases, these 
quasi-lD regions grow both in width and length and their 
correlations ultimately cross over to a pure 2D character. 
This dimensionality effect must be taken into account 
when using finite size systems to infer the critical behav- 
ior. 

In summary, we have addressed finite temperature 
properties of inhomogeneous Fermi-Hubbard systems us- 
ing an ab-initio approach. Our findings of enhanced an- 
tiferromagnetic and pairing correlations just below the 
temperature scale T ^ t thus open important perspec- 
tives for current experiments with ultracold fermions in 
optical lattices. While the lowest temperatures reported 
here are well above the c?-wave coherence temperature, 
the enhanced signal in local properties is a promising 
sign: these results and the purity of the experimental 
optical lattice suggest that, in contrast to solid state sys- 
tems, temperatures of the order of the hopping scale may 
be enough to observe clear local signature of both mag- 
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Shortest non-trivial correlations If it would be possi- 
ble to address low enough temperatures that the various 
correlations in the trap would all saturate to their largest 
value, we would then use an order parameter to establish 
where and how the system is most likely to order. Exper- 
iments are far from that regime and we therefore focus 
on the task on selecting short-range correlations which 
are representative (in the sense of having the right spa- 
tial, temperature and U dependence) of the behavior of 
a true order parameter. 

A useful analogy is provided by considering the min- 
imal length scale at which the spin correlation function 
reflects the exchange energy J = At^ /U . Clearly the 
temperature dependence of the shortest separation spin 
correlator. Coo (the local moment) is completely con- 
trolled by the energy scale U rather than J, as can 
be seen by rewriting it in terms of density properties. 
Coo = (") - 2 

To develop the same sort of analysis in the case of 
pairing correlations, consider now the contribution that 
comes from creating and destroying a pair on the same 




bond. One can derive the following identity, 

(A,, At^.), = - (S, . S,) + l(Ar.,. An,) - 1 ^ {c. J^^ ■ 

a 

(1) 

The three contributions to (AijA|j)c are reported in 
Fig. 1 as a function of chemical potential for a homo- 
geneous 8 X 8-site cluster, U = 6t and T = 0.5 1. One 
can see that the largest contribution is given by Si • Sj, 
especially in the most interesting density regimes at and 
near half-filling. This simple analysis reveals that not 
only is the energy scale trivially inherited from that of 
spin fluctuations but also that the location of the maxi- 
mum is pinned where • Sj is bigger, i.e. at half-filling. 

When using the definition of the spatially dependent 
d-wave pairing function, the local and nearest-neighbor 
contributions can be graphically pictured as in panels (a) 
and (b) of Fig. 2, and they therefore respectively receive 
four and one contributions of the pathological type writ- 
ten in Eq. (1). These pairing profiles are given in panels 
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FIG. 1: The three contributions to the connected pairing cor- 
relation function that corresponds to creating and destroying 
a pair of particles on the same (nearest-neighbor) bond as 
a function of chemical potential for a homogeneous system 
with U = 6t and T = 0.5 1. fi — corresponds to half-filling. 
The large contribution from (Si ■ Sj) suggests that this con- 
tribution will dominate the shortest range pair correlations 
and hence that the dependence on temperature and [/ in a 
trapped system of somewhat longer range Pxyii) should be 
studied to get the clearest insight into incipient pair order. 



FIG. 2: Panels (a) and (b) report a schematic representation 
of the local and nearest-neighbor pair overlap defining Poo(i) 
and Pio(i). DQMC values for these quantities are plotted in 
panels (c) and (d). Their shape as a function of position is 
very similar and strongly reminiscent of Cio(i). The drop in 
amplitude perfectly correlates to the decrease in the number 
of overlapping bonds in panels (a) and (b). As previously 
suggested by Fig. 1, both quantities are strongly dominated 
by Cio(i) making it difficult to make use of them as a measure 
of the strength of pairing correlations. 
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FIG. 3: Solid symbols: pairing amplitude in the d (or- 
ange/triangle) and extended-s (green/circle) channels at 
T = 0.31*. Empty symbols: pairing amplitude in the d 
(red/triangle) and extended-s (blue/circle) channels at T = 
0.54t. Different symmetry of the order parameter are peaked 
in different regions. 

(c) and (d) of Fig. 2: their shape mimics the Cio(i) profile 
and the intensity decreases by a factor of four, reflecting 
the ratio of overlapping bonds in the upper panels. It 
is thus clear that such correlations cannot be used as a 
meaningful measure of pairing in the system and that the 
shortest non-trivial correlations are Cio(«) and Pii(i). 

Extended s-wave pairing The extended s-wave pair 
creation operator is defined by replacing in the main 
text by: 

Af = ^(Al,,+, + At + A[,+ - + A[,_.). (2) 

Fig. 3 reports the spatial variation of the (1,1) cor- 
relation for both extended- s-wave and d-wave type at 
T = 0.54* and T = OMt. While both correlations in- 
crease with decreasing temperature in regions surround- 
ing the Mott insulating plateaux (whose radius is about 8 
lattice spacings) the d-wave is peaked in narrow regions 
centered at density n = 1.0 ± 0.2 while the s-wave in 
broader ones around n = 1.0 ± 0.4. This is in agreement 
with existing ground state calculations on related homo- 
geneous models in small clusters where the extended s- 



wave order parameter was found to be larger than the 
one of d-wave order in the overdoped regime [1, 2]. Our 
present finite temperature analysis reveals that if an in- 
stability were to occur, it would start developing around 
n = 0.8 for d-wave pairing and around n = 0.6 for ex- 
tended s-wave pairing. 

Calculation of the entropy The free energy of a non- 
interacting system can be computed using the identity 

F{a, T) - F{Q, T) = ^^^^d^' (3) 

where F{X, T) is the free energy associated to the hamil- 
tonian H{X) = Hq + XHi. 

For the Hubbard model we can set X ~ U and Hi = 
ni^Uii so that the free energy can be expressed as an 
integral of the double occupancy 

F{U,T)-F{Q,T)^ dU y^{n,^n,^)u.T- (4) 
Jo 

In this paper we have used the LDA to estimate 
the value of the local double occupancies since ab-initio 
DQMC results at the same temperature have confirmed 
that such a quantity can be reliably described in this 
approximation. 

We determine the entropy using the thermodynamic 
identity S = /3{E — F) with E, once again, obtained 
within LDA. Such an estimate of E is certainly of in- 
ferior accuracy if compared to (n-fnj,), due to stronger 
sensitivity of the kinetic energy on boundary conditions. 
However, for those U values for which we have carried out 
full DQMC simulations in a trap, we found the accuracy 
in estimating E to be in the 5% range and completely 
adequate for the kind of conclusions we are drawing in 
the main body of this work. 
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